function output = GeneratePO(~,varargin)
options_fmincon = optimoptions('fmincon', 'MaxFunctionEvaluations', 20000, ...
    'MaxIterations', 3000,'ConstraintTolerance',5e-5,'Algorithm','Interior-point','stepTolerance',1e-21);
options=odeset('RelTol',1e-13,'AbsTol',1e-13);

%% Input List
% 1. System: EM (Earth-Moon), SE (Sun-Earth)
%
% 2. PO_Location: {L1,L2,L4,L5}
%     Location of periodic orbit of interest
%
% 3. PO_Type: {Lyapunov,Halo_Sourn,Halo_Northern,Verthetical}
%     Possible Inputs:
%     L1&2: Lyapunov,Halo_Southern,Halo_Northern
%     L4&5: Lyapunov,Vertical
%
% 4. JacobiTarget: Target Jacobi's integral
%     J = single floating Number

%%
numvarargs = length(varargin);
% Number of optional arguments

if mod(numvarargs, 2) == 1
    error('Error: Please check name-value pair arguments');
end

name_arguments = varargin(1:2:end);
value_arguments = varargin(2:2:end);
% Initialize name-value arguments

for ii = 1:length(name_arguments)
    switch lower(name_arguments{ii})
        case 'massratio'
            u = value_arguments{ii};
        case 'system'
            System_Name = value_arguments{ii};
        case 'polocation'
            PO_location = value_arguments{ii};
        case 'potype'
            PO_Type = value_arguments{ii};
        case 'jacobitarget'
            J_Target = value_arguments{ii};
    end
end

FileName = append(System_Name,PO_location,'_',PO_Type,'.mat');
load(FileName);

%% Determine Feasibility

N = length(IC);
%[] Number of initial conditions

J_list = zeros(1,N);
%[] Pre-allocate Jacobi's Integral value

for ii = 1:N

    J_list(ii) = jacobiConst(IC(1:3,ii),IC(4:6,ii),u);
    % Determine Jacobi's integral of all initial conditions
end

J_min = min(J_list);
J_max = max(J_list);
fprintf('Min J: %f\n',J_min)
fprintf('Max J: %f\n',J_max)
% Determine min and maximum jacobi's integral in the list

if or(J_Target<J_min,J_Target>J_max);

    error('JACOBI INTEGRAL TARGET IS OUT OF MIN/MAX BOUNDARY OF INITIAL CONDITION MAT FILE\n Minimum possible J=%3.5f\n Maximum possible J=%3.5f\n\n',J_min,J_max')
end
% Catch whether jacobi's integral is within the region of provided initial
% guesses

%% Select PO closest to J_target
J_temp = abs(J_list-J_Target);
ind_Closest_PO = find(J_temp == min(J_temp));
%[] Index of best PO

So = IC(:,ind_Closest_PO);
FT = T(ind_Closest_PO);

%% Fine PO of interest
PO_Location_Type = append(PO_location,PO_Type);
fprintf('Runing FMINCON\n\n')
switch lower(PO_Location_Type)
    case 'l1lyapunov'
        % Design variables = [x+dx;0;0;0;vy+dvy;0]
        Xo = [So(1),So(5),FT];
        [Xopt] = fmincon(@CostFunction,...
        Xo,...     % Initial guess
        [], [],... % Linear inequality constraint
        [], [],... % Linear equality constraint
        [], [],... % Boundary constraint
        @(X)DynamicalConst_Lyapunov(X,u,options,J_Target),... % Nonlinear Constraint and options.
        options_fmincon);   % Options for fmincon
        output.X = [Xopt(1);0;0;0;Xopt(2);0];
        output.FT = 2*Xopt(3);
    case 'l1northern'
        % Design variables = [x+dx;0;z+dz;0;vy+dvy;0]
        Xo = [So(1),So(3),So(5),FT];
        [Xopt] = fmincon(@CostFunction,...
        Xo,...     % Initial guess
        [], [],... % Linear inequality constraint
        [], [],... % Linear equality constraint
        [], [],... % Boundary constraint
        @(X)DynamicalConst_Halo(X,u,options,J_Target),... % Nonlinear Constraint and options.
        options_fmincon);   % Options for fmincon
        output.X = [Xopt(1);0;Xopt(2);0;Xopt(3);0];
        output.FT = 2*Xopt(4);
    case 'l1southern'
        % Design variables = [x+dx;0;z+dz;0;vy+dvy;0]
        Xo = [So(1),So(3),So(5),FT];
        [Xopt] = fmincon(@CostFunction,...
        Xo,...     % Initial guess
        [], [],... % Linear inequality constraint
        [], [],... % Linear equality constraint
        [], [],... % Boundary constraint
        @(X)DynamicalConst_Halo(X,u,options,J_Target),... % Nonlinear Constraint and options.
        options_fmincon);   % Options for fmincon
        output.X = [Xopt(1);0;Xopt(2);0;Xopt(3);0];
        output.FT = 2*Xopt(4);

    case 'l2lyapunov'
        % Design variables = [x+dx;0;0;0;vy+dvy;0]
        Xo = [So(1),So(5),FT];
        [Xopt] = fmincon(@CostFunction,...
        Xo,...     % Initial guess
        [], [],... % Linear inequality constraint
        [], [],... % Linear equality constraint
        [], [],... % Boundary constraint
        @(X)DynamicalConst_Lyapunov(X,u,options,J_Target),... % Nonlinear Constraint and options.
        options_fmincon);   % Options for fmincon
        output.X = [Xopt(1);0;0;0;Xopt(2);0];
        output.FT = 2*Xopt(3);
    case 'l2northern'
        % Design variables = [x+dx;0;z+dz;0;vy+dvy;0]
        Xo = [So(1),So(3),So(5),FT];
        [Xopt] = fmincon(@CostFunction,...
        Xo,...     % Initial guess
        [], [],... % Linear inequality constraint
        [], [],... % Linear equality constraint
        [], [],... % Boundary constraint
        @(X)DynamicalConst_Halo(X,u,options,J_Target),... % Nonlinear Constraint and options.
        options_fmincon);   % Options for fmincon
        output.X = [Xopt(1);0;Xopt(2);0;Xopt(3);0];
        output.FT = 2*Xopt(4);
    case 'l2southern'
        % Design variables = [x+dx;0;z+dz;0;vy+dvy;0]
        Xo = [So(1),So(3),So(5),FT];
        [Xopt] = fmincon(@CostFunction,...
        Xo,...     % Initial guess
        [], [],... % Linear inequality constraint
        [], [],... % Linear equality constraint
        [], [],... % Boundary constraint
        @(X)DynamicalConst_Halo(X,u,options,J_Target),... % Nonlinear Constraint and options.
        options_fmincon);   % Options for fmincon
        output.X = [Xopt(1);0;Xopt(2);0;Xopt(3);0];
        output.FT = 2*Xopt(4);
    case 'l4lyapunov'
        % Design variables = [x+dx;0.866025403784439;0;vx+dvx;vy+dvy;0]
        Xo = [So(1),So(4),So(5),FT];
    case 'l5lyapunov'
        Xo = [So(1),So(4),So(5),FT];
end

end

function J = CostFunction(X)
J = 0;
end


function [c,ceq] = DynamicalConst_Lyapunov(X,u,options,J_target)
% Design variables = [x+dx;0;0;0;vy+dvy;0]
So = [X(1);0;0;0;X(2);0];
FT = [0,2*X(3)];
[~,S] = ode113(@(t,S)CR3BP_n(t,S,u),FT,So,options); S = S';
% plot3(S(1,:),S(2,:),S(3,:))
Jacobi = jacobiConst(So(1:3),So(4:6),u);
c = [];
ceq = [S(1:3,end)-So(1:3);...
       S(4:6,end)-So(4:6);...
       Jacobi-J_target];
end

function [c,ceq] = DynamicalConst_Halo(X,u,options,J_target)
% Design variables = [x+dx;0;z+dz;0;vy+dvy;0]
So = [X(1);0;X(2);0;X(3);0];
FT = [0,2*X(4)];
[~,S] = ode113(@(t,S)CR3BP_n(t,S,u),FT,So,options); S = S';
% plot3(S(1,:),S(2,:),S(3,:))
Jacobi = jacobiConst(So(1:3),So(4:6),u);
c = [];
ceq = [S(1:3,end)-So(1:3);...
       S(4:6,end)-So(4:6);...
       Jacobi-J_target];
end

% function [c,ceq] = DynamicalConst_L4Lyapunov(X,u,options,J_target)
% end
% 
% function [c,ceq] = DynamicalConst_L5Lyapunov(X,u,options,J_target)
% end